#include "T2KParticle.h"

class NeutrinoUtils
{
public:
    enum NeutrinoParent{
        kNull = -999,
        kBad = -2,
        kLost = -1,
        kOther=0,
        kPion,
        kKaon,
        kMuon
    };
    
    static void IniCC();
    static void IniCCQE();
    static void IniCC1P1Pi();
    
    //--- Reconstruction related
    static const TLorentzVector * fNeutrinoRec;
    static T2KParticle * fMuonRec;
    static T2KParticle * fProtonRec;
    static T2KParticle * fPionRec;
    static T2KParticle * fMuonGuess;
    static const TLorentzVector * fNeutrinoGuess;
    
    //--- Simulation related
    static T2KParticle * fMuonSim;
    static T2KParticle * fProtonSim;
    static T2KParticle * fPionSim;
    static T2KParticle * fNeutronToy;
    static T2KParticle * fProtonToy;
    static T2KParticle * fProtonDet;
    
    static const TLorentzVector * fNeutrinoSim;
    
    //--- Transverse related
    static TVector3 * fDeltaRecPt;
    static Double_t * fpTTRec;
    static Double_t * fAbspTTRec;
    static Double_t * fpTTSim;
    static Double_t * fAbspTTSim;
    static Double_t * fDeltapTTRec;
    static Double_t * fAbsDeltapTTRec;
    static Double_t * fDeltapTTSim;
    static Double_t * fAbsDeltapTTSim;
    
    static TVector3 * fDeltaSimPt;
    static TVector3 * fDeltaToyPt;
    static TVector3 * fDeltaDetPt;
    
    
    
    static const TVector3 * fNeutrinoParentDecPointSim;
    
    static Double_t fOsciMuMuProb;
    static Int_t fkOsciMuSurvive;
    static Int_t fNeutrinoParentPDG;
    static Int_t fNeutrinoParentType;
    static Int_t fNeutrinoType;
    static Int_t fNeutMode;
    
    static void SetMuonNRec(const Int_t var){ fMuonNRec = var; }
    static Int_t GetMuonNRec(){return fMuonNRec;}
    static void SetMuonNSim(const Int_t var){ fMuonNSim = var; }
    static Int_t GetMuonNSim(){return fMuonNSim;}
    
    static void SetProtonNRec(const Int_t var){ fProtonNRec = var; }
    static Int_t GetProtonNRec(){return fProtonNRec;}
    static void SetProtonNSim(const Int_t var){ fProtonNSim = var; }
    static Int_t GetProtonNSim(){return fProtonNSim;}
    
    static void SetMultiplicity(const Int_t mul){fMultiplicity = mul;}
    static Int_t GetMultiplicity(){return fMultiplicity;}
    
    //only rec. meaningful; phi angle in the transver plane w.r.t. beam direction, -180 ~ 180, by definition muon at (-)180.
    static Double_t fBetaRec[10];
    static Double_t fGammaRec[10];
    
    //-------- Calculations ---------
    static void SetDeltaPt();
    //PTT VARIABLE
    static void SetpTT();
    
    //--- Reconstruction related
    Double_t CalcBetaRec(const Double_t dir0, const Double_t dir1, const Double_t dir2);
    Double_t CalcGammaRec(const Double_t dir0, const Double_t dir1, const Double_t dir2);
    
    static void SetRecGuessKinematics(const Bool_t kdummy=kFALSE);
    
    //---Simulation related
    static void MuonNSimSanityCheck();
    
    static void SetSimKinematics();
    
    //-------- IO ---------
    static TTree * GetTree();
    static TList * GetHistList();
    
    static Int_t FillCount(TList *lin, const Int_t var){ return NeutrinoTools::FillTH1I(lin, kcount, var);}
    
private:
    //==========================================================
    //Variables and setters
    //==========================================================
    //--- Reconstruction related
    static Double_t fdNeutrinoGuessNeutrinoSim;
    static Double_t fdNeutrinoGuessMuonRec;
    static const Int_t fCutNR;
    static Double_t fGuessCorrection;
    static void SetdNeutrinoGuessNeutrinoSim();
    static void SetdNeutrinoGuessMuonRec();
    
    static Int_t fMuonNRec;
    static Int_t fProtonNRec;
    
    static Int_t fMultiplicity;
    static Double_t fMeanBetaRec;
    static void SetMeanBetaRec(const Bool_t kdummy = kFALSE);
    static Double_t fMeanGammaRec;//use Mean instead of sum to get rid of multiplicity
    static void SetMeanGammaRec(const Bool_t kdummy = kFALSE);
    
    ///---Simulation related
    static Bool_t IsCC(){ return TMath::Abs(fNeutMode)<30;}
    
    static Int_t fMuonNSim;
    static Int_t fProtonNSim;
    
    enum histID{
        kcount,
        kDUMMY
    };
    
};

//__________________________________________________________

const TLorentzVector * NeutrinoUtils::fNeutrinoRec = 0x0;
T2KParticle *NeutrinoUtils::fMuonRec=0x0;
T2KParticle *NeutrinoUtils::fProtonRec=0x0;
T2KParticle *NeutrinoUtils::fPionRec=0x0;
TVector3 * NeutrinoUtils::fDeltaRecPt = 0x0;

//Transverse related
Double_t * NeutrinoUtils::fpTTRec = 0x0;
Double_t * NeutrinoUtils::fAbspTTRec = 0x0;
Double_t * NeutrinoUtils::fDeltapTTRec = 0x0;
Double_t * NeutrinoUtils::fAbsDeltapTTRec = 0x0;
Double_t * NeutrinoUtils::fDeltapTTSim = 0x0;
Double_t * NeutrinoUtils::fAbsDeltapTTSim = 0x0;

const TLorentzVector * NeutrinoUtils::fNeutrinoGuess = 0x0;
T2KParticle *NeutrinoUtils::fMuonGuess=0x0;
const TLorentzVector * NeutrinoUtils::fNeutrinoSim = 0x0;
T2KParticle *NeutrinoUtils::fMuonSim=0x0;
T2KParticle *NeutrinoUtils::fProtonSim=0x0;
T2KParticle *NeutrinoUtils::fPionSim=0x0;
T2KParticle *NeutrinoUtils::fNeutronToy=0x0;
T2KParticle *NeutrinoUtils::fProtonToy=0x0;
T2KParticle *NeutrinoUtils::fProtonDet=0x0;
TVector3 * NeutrinoUtils::fDeltaSimPt = 0x0;
TVector3 * NeutrinoUtils::fDeltaToyPt = 0x0;
TVector3 * NeutrinoUtils::fDeltaDetPt = 0x0;
Double_t * NeutrinoUtils::fpTTSim = 0x0;
Double_t * NeutrinoUtils::fAbspTTSim = 0x0;
const TVector3 * NeutrinoUtils::fNeutrinoParentDecPointSim = 0x0;

Double_t NeutrinoUtils::fOsciMuMuProb = NeutrinoTools::kIniValue;
Int_t NeutrinoUtils::fkOsciMuSurvive = NeutrinoTools::kIniValue;
Int_t NeutrinoUtils::fNeutrinoParentPDG = NeutrinoUtils::kNull;
Int_t NeutrinoUtils::fNeutrinoParentType = NeutrinoUtils::kNull;
Int_t NeutrinoUtils::fNeutrinoType = NeutrinoTools::kIniValue;
Int_t NeutrinoUtils::fNeutMode = NeutrinoTools::kIniValue;

Int_t NeutrinoUtils::fMuonNRec = NeutrinoTools::kIniValue;
Int_t NeutrinoUtils::fMuonNSim = NeutrinoTools::kIniValue;
Int_t NeutrinoUtils::fProtonNRec = NeutrinoTools::kIniValue;
Int_t NeutrinoUtils::fProtonNSim = NeutrinoTools::kIniValue;

Int_t NeutrinoUtils::fMultiplicity = NeutrinoTools::kIniValue;

Double_t NeutrinoUtils::fBetaRec[10]={NeutrinoTools::kIniValue, NeutrinoTools::kIniValue, NeutrinoTools::kIniValue, NeutrinoTools::kIniValue, NeutrinoTools::kIniValue, NeutrinoTools::kIniValue,NeutrinoTools::kIniValue, NeutrinoTools::kIniValue, NeutrinoTools::kIniValue,NeutrinoTools::kIniValue};
Double_t NeutrinoUtils::fMeanBetaRec = NeutrinoTools::kIniValue;
Double_t NeutrinoUtils::fGammaRec[10]={NeutrinoTools::kIniValue, NeutrinoTools::kIniValue, NeutrinoTools::kIniValue, NeutrinoTools::kIniValue, NeutrinoTools::kIniValue, NeutrinoTools::kIniValue,NeutrinoTools::kIniValue, NeutrinoTools::kIniValue, NeutrinoTools::kIniValue,NeutrinoTools::kIniValue};
Double_t NeutrinoUtils::fMeanGammaRec = NeutrinoTools::kIniValue;

Double_t NeutrinoUtils::fdNeutrinoGuessNeutrinoSim = NeutrinoTools::kIniValue;
Double_t NeutrinoUtils::fdNeutrinoGuessMuonRec = NeutrinoTools::kIniValue;

Double_t NeutrinoUtils::fGuessCorrection = 1;
const Int_t NeutrinoUtils::fCutNR = 1;

//_______________________________________________________________________________________________
//_______________________________________________________________________________________________
void NeutrinoUtils::IniCC()
{
    delete NeutrinoUtils::fMuonSim; NeutrinoUtils::fMuonSim= new T2KParticle;
    delete NeutrinoUtils::fMuonRec; NeutrinoUtils::fMuonRec= new T2KParticle;
    delete NeutrinoUtils::fMuonGuess; NeutrinoUtils::fMuonGuess= new T2KParticle;
}

void NeutrinoUtils::IniCCQE()
{
    IniCC();
    delete NeutrinoUtils::fProtonSim; NeutrinoUtils::fProtonSim= new T2KParticle;
    delete NeutrinoUtils::fNeutronToy; NeutrinoUtils::fNeutronToy= new T2KParticle;
    delete NeutrinoUtils::fProtonToy; NeutrinoUtils::fProtonToy= new T2KParticle;
    delete NeutrinoUtils::fProtonDet; NeutrinoUtils::fProtonDet= new T2KParticle;
    delete NeutrinoUtils::fProtonRec; NeutrinoUtils::fProtonRec= new T2KParticle;
    delete NeutrinoUtils::fDeltaRecPt; NeutrinoUtils::fDeltaRecPt=new TVector3;
    delete NeutrinoUtils::fDeltaSimPt; NeutrinoUtils::fDeltaSimPt=new TVector3;
    delete NeutrinoUtils::fDeltaToyPt; NeutrinoUtils::fDeltaToyPt=new TVector3;
    delete NeutrinoUtils::fDeltaDetPt; NeutrinoUtils::fDeltaDetPt=new TVector3;
    //Additional code for double transverse
    delete NeutrinoUtils::fpTTRec; NeutrinoUtils::fpTTRec=new Double_t;
    delete NeutrinoUtils::fpTTSim; NeutrinoUtils::fpTTSim=new Double_t;
    delete NeutrinoUtils::fAbspTTRec; NeutrinoUtils::fAbspTTRec=new Double_t;
    delete NeutrinoUtils::fAbspTTSim; NeutrinoUtils::fAbspTTSim=new Double_t;
}

void NeutrinoUtils::IniCC1P1Pi()
{
    IniCCQE();
    delete NeutrinoUtils::fPionSim; NeutrinoUtils::fPionSim= new T2KParticle;
    delete NeutrinoUtils::fPionRec; NeutrinoUtils::fPionRec= new T2KParticle;
    delete NeutrinoUtils::fDeltapTTRec; NeutrinoUtils::fDeltapTTRec= new T2KParticle;
    delete NeutrinoUtils::fAbsDeltapTTRec; NeutrinoUtils::fAbsDeltapTTRec= new T2KParticle;
    delete NeutrinoUtils::fDeltapTTSim; NeutrinoUtils::fDeltapTTSim= new T2KParticle;
    delete NeutrinoUtils::fAbsDeltapTTSim; NeutrinoUtils::fAbsDeltapTTSim=new T2KParticle;
}

void NeutrinoUtils::SetRecGuessKinematics(const Bool_t kdummy)
{
    SetdNeutrinoGuessMuonRec();
    SetMeanBetaRec(kdummy);
    SetMeanGammaRec(kdummy);
}

void NeutrinoUtils::SetSimKinematics()
{
    SetdNeutrinoGuessNeutrinoSim();
}

void NeutrinoUtils::SetDeltaPt()
{
    delete fDeltaRecPt; fDeltaRecPt = new TVector3;
    delete fDeltaSimPt; fDeltaSimPt = new TVector3;
    delete fDeltaToyPt; fDeltaToyPt = new TVector3;
    delete fDeltaDetPt; fDeltaDetPt = new TVector3;
    
    if(fMuonRec->fPt && fProtonRec->fPt){
        //has to be ordered: muon, then proton
        NeutrinoTools::SetDeltaPt(fDeltaRecPt, fMuonRec->fPt, fProtonRec->fPt);
    }
    if(fMuonSim->fPt && fProtonSim->fPt){
        NeutrinoTools::SetDeltaPt(fDeltaSimPt, fMuonSim->fPt, fProtonSim->fPt);
    }
    
    if(fMuonSim->fPt && fProtonToy->fPt){
        NeutrinoTools::SetDeltaPt(fDeltaToyPt,  fMuonSim->fPt, fProtonToy->fPt);
    }
    
    if(fMuonSim->fPt && fProtonDet->fPt){
        NeutrinoTools::SetDeltaPt(fDeltaDetPt,  fMuonSim->fPt, fProtonDet->fPt);
    }
}

void NeutrinoUtils::SetpTT()
{
    delete fpTTRec; fpTTRec = new Double_t;
    delete fpTTSim; fpTTSim = new Double_t;
    
    if(fNeutrinoRec && fMuonRec->fPt && fProtonRec->fPt){
        NeutrinoTools::SetpTT(fpTTRec, fNeutrinoRec, fMuonRec->fPt, fProtonRec->fPt);
        NeutrinoTools::SetAbspTT(fAbspTTRec, fNeutrinoRec, fMuonRec->fPt, fProtonRec->fPt);
    }
    if(fNeutrinoSim && fMuonSim->fPt && fProtonSim->fPt){
        NeutrinoTools::SetpTT(fpTTSim, fNeutrinoSim, fMuonSim->fPt, fProtonSim->fPt);
        NeutrinoTools::SetAbspTT(fAbspTTSim, fNeutrinoSim, fMuonSim->fPt, fProtonSim->fPt);
        //NeutrinoTools::SetpTT(fpTTSim, fNeutrinoSim, fMuonSim->fPt, fProtonDet->fPt);
        //NeutrinoTools::SetAbspTT(fAbspTTSim, fNeutrinoSim, fMuonSim->fPt, fProtonDet->fPt);
    }
}

void NeutrinoUtils::SetDeltapTT(){
    delete fDeltapTTRec; fDeltapTTRec = new Double_t;
    delete fDeltapTTSim; fDeltapTTSim = new Double_t;
    delete fAbsDeltapTTRec; fAbsDeltapTTRec = new Double_t;
    delete fAbsDeltapTTSim; fAbsDeltapTTSim = new Double_t;
    
    if(fNeutrinoRec && fMuonRec->fPt && fProtonRec->fPt && fPionRec->fPt){
        NeutrinoTools::SetDeltapTT(fDeltapTTRec, fNeutrinoRec, fMuonRec->fPt, fProtonRec->fPt, fPionRec->fPt);
        NeutrinoTools::SetDeltaAbspTT(fAbsDeltapTTRec, fNeutrinoRec, fMuonRec->fPt, fProtonRec->fPt, fPionRec->fPt);
        
    }
    if(fNeutrinoSim && fMuonSim->fPt && fProtonSim->fPt && fPionSim->fPt){
        NeutrinoTools::SetDeltapTT(fDeltapTTSim, fNeutrinoSim, fMuonSim->fPt, fProtonSim->fPt, fPionSim->fPt);
        NeutrinoTools::SetDeltaAbspTT(fAbsDeltapTTSim, fNeutrinoSim, fMuonSim->fPt, fProtonSim->fPt, fPionSim->fPt);
    }
}

void NeutrinoUtils::SetdNeutrinoGuessNeutrinoSim()
{
    //flat in cos, not tha angle
    //tested with = fNeutrinoRec; should be 0; confirmed
    if(!fNeutrinoGuess){
        fdNeutrinoGuessNeutrinoSim = NeutrinoTools::kNullPointer;
    }
    else{
        const Double_t tmpcos= NeutrinoTools::GetCos(fNeutrinoGuess, fNeutrinoSim, "NeutrinoUtils::SetdNeutrinoGuessNeutrinoSimGetCos");
        if(tmpcos== NeutrinoTools::kZeroDivider){
            fdNeutrinoGuessNeutrinoSim = NeutrinoTools::kZeroDivider;
        }
        else{
            fdNeutrinoGuessNeutrinoSim =  tmpcos - 1;
        }
    }
}

void NeutrinoUtils::SetdNeutrinoGuessMuonRec()
{
    //flat in cos, not tha angle
    
    //calculation confirmed by
    //tree->Draw("fdNeutrinoGuessMuonRec","fMultiplicity==1") //should = 1
    
    //tested with = &(fMuonRec->Vect())// should be 0; confirmed
    if(!fMuonRec->fLorentzVec || !fNeutrinoGuess){
        fdNeutrinoGuessMuonRec = NeutrinoTools::kNullPointer;
    }
    else{
        const Double_t tmpcos = NeutrinoTools::GetCos(fNeutrinoGuess, fMuonRec->fLorentzVec, "NeutrinoUtils::SetdNeutrinoGuessMuonRecGetCos");
        if(tmpcos== NeutrinoTools::kZeroDivider){
            fdNeutrinoGuessMuonRec = NeutrinoTools::kZeroDivider;
        }
        else{
            fdNeutrinoGuessMuonRec = tmpcos - 1;
        }
    }
}

void NeutrinoUtils::SetMeanBetaRec(const Bool_t kdummy)
{
    if(kdummy){
        fMeanBetaRec = NeutrinoTools::kIniValue;
        return;
    }
    
    fMeanBetaRec = 0;
    for(Int_t ipar=0; ipar<fMultiplicity; ipar++){
        //remove muon contribution
        if(ipar == fMuonRec->fID){
            continue;
        }
        
        fMeanBetaRec += fBetaRec[ipar];
    }
    
    if(fMultiplicity==1){
        fMeanBetaRec=NeutrinoTools::kZeroDivider;
    }
    else{//multiplicity >= 2
        fMeanBetaRec /= (fMultiplicity-1);
    }
}

void NeutrinoUtils::SetMeanGammaRec(const Bool_t kdummy)
{
    if(kdummy){
        fMeanGammaRec = NeutrinoTools::kIniValue;
        return;
    }
    
    fMeanGammaRec = 0;
    for(Int_t ipar=0; ipar<fMultiplicity; ipar++){
        //remove muon contribution
        if(ipar == fMuonRec->fID){
            continue;
        }
        fMeanGammaRec += fGammaRec[ipar];
    }
    
    if(fMultiplicity==1){
        fMeanGammaRec = NeutrinoTools::kZeroDivider;
    }
    else{
        fMeanGammaRec /= (fMultiplicity-1);
    }
}

void NeutrinoUtils::MuonNSimSanityCheck()
{
    //sanity check -1
    if(fMuonNSim==0 && IsCC()){
        printf("MuonNSimSanityCheck wrong fMuonNSim CC %d %d\n", fMuonNSim, fNeutMode); exit(1);
    }
    //sanity check -2
    if(fMuonNSim>0 && !IsCC()){
        printf("MuonNSimSanityCheck wrong fMuonNSim NC %d %d\n", fMuonNSim, fNeutMode); exit(1);
    }
}

Double_t NeutrinoUtils::CalcBetaRec(const Double_t dir0, const Double_t dir1, const Double_t dir2)
{
    //
    //magnitude does not magger
    //w.r.t. beam direction (z-axis) and muon opposite direction  (x-axis)
    //-180 ~ 180
    //
    if(!fMuonRec->fPt){
        return NeutrinoTools::kNullPointer;
    }
    
    const TVector3 *vt = NeutrinoTools::GetVecT(fNeutrinoRec, dir0, dir1, dir2);
    
    //muon pt direction := -180 deg, i.e. -x axis
    const TVector3 xaxis = (*fMuonRec->fPt)*(-1);
    
    //0 ~ 180
    
    
    Double_t beta = NeutrinoTools::GetAngle(&xaxis, vt, "NeutrinoUtils::CalcBetaRec");
    
    if( (fNeutrinoRec->Vect()).Dot(xaxis.Cross(*vt)) < 0 )
        beta *= -1;
    
    if(beta<-180 + 1e-3){
        beta+= 360;
    }
    
    delete vt;
    return beta;
}

Double_t NeutrinoUtils::CalcGammaRec(const Double_t dir0, const Double_t dir1, const Double_t dir2)
{
    const TVector3 vec(dir0, dir1, dir2);
    
    const TVector3 refdir=fNeutrinoRec->Vect();
    return NeutrinoTools::GetAngle(&vec, &refdir, "NeutrinoUtils::CalcGammaRec");
}

//___________________________________________________________________________________

//___________________________________________________________________________________________________

TTree * NeutrinoUtils::GetTree()
{
    TTree * tout = new TTree("tree","tree");
    
    const Int_t spl = 1;
    
    fMuonRec->SetTree("fMuonRec",tout,spl);
    fMuonSim->SetTree("fMuonSim",tout,spl);
    
    if(fProtonRec && fProtonSim){
        fProtonRec->SetTree("fProtonRec",tout,spl);
        fProtonSim->SetTree("fProtonSim",tout,spl);
    }
    
    if(fDeltaRecPt && fDeltaSimPt){
        tout->Branch("fDeltaRecPt", "TVector3", &(fDeltaRecPt), 128000, spl);
        tout->Branch("fDeltaSimPt", "TVector3", &(fDeltaSimPt), 128000, spl);
    }
    
    if(fpTTRec && fpTTSim){
        tout->Branch("fpTTRec",&(*fpTTRec),"fpTTRec/D");
        tout->Branch("fpTTSim",&(*fpTTSim),"fpTTSim/D");
    }
    
    if(fAbspTTRec && fAbspTTSim){
        tout->Branch("fAbspTTRec",&(*fAbspTTRec),"fAbspTTRec/D");
        tout->Branch("fAbspTTSim",&(*fAbspTTSim),"fAbspTTSim/D");
    }
    
    if(fDeltapTTRec && fDeltapTTSim){
        tout->Branch("fDeltapTTRec",&(*fDeltapTTRec),"fDeltapTTRec/D");
        tout->Branch("fDeltapTTSim",&(*fDeltapTTSim),"fDetlapTTSim/D");
    }
    if(fAbsDeltapTTRec && fAbsDeltapTTSim){
        tout->Branch("fAbsDeltapTTRec",&(*fAbsDeltapTTRec),"fAbsDeltapTTRec/D");
        tout->Branch("fAbsDeltapTTSim",&(*fAbsDeltapTTSim),"fAbsDeltapTTSim/D");
    }
    
    if(fProtonToy && fDeltaToyPt){
        fNeutronToy->SetTree("fNeutronToy",tout,spl);
        fProtonToy->SetTree("fProtonToy",tout,spl);
        tout->Branch("fDeltaToyPt", "TVector3", &(fDeltaToyPt), 128000, spl);
    }
    
    if(fProtonDet && fDeltaDetPt){
        fProtonDet->SetTree("fProtonDet",tout,spl);
        tout->Branch("fDeltaDetPt", "TVector3", &(fDeltaDetPt), 128000, spl);
    }
    
    tout->Branch("fNeutrinoSim","TLorentzVector", &fNeutrinoSim, 128000, spl);
    tout->Branch("fNeutrinoRec", "TVector3", &(fMuonRec->fNuDir), 128000, spl);
    tout->Branch("fNeutrinoGuess", "TVector3", &(fMuonGuess->fNuDir), 128000, spl);
    
    tout->Branch("fdNeutrinoGuessNeutrinoSim", &fdNeutrinoGuessNeutrinoSim, "fdNeutrinoGuessNeutrinoSim/D");
    tout->Branch("fdNeutrinoGuessMuonRec", &fdNeutrinoGuessMuonRec, "fdNeutrinoGuessMuonRec/D");
    
    tout->Branch("fGuessCorrection", &fGuessCorrection, "fGuessCorrection/D");
    tout->Branch("fMuonPtGuess", &(fMuonGuess->fPtMag), "fMuonPtGuess/D");
    
    tout->Branch("fOsciMuMuProb",&fOsciMuMuProb,"fOsciMuMuProb/D");
    tout->Branch("fkOsciMuSurvive",&fkOsciMuSurvive,"fkOsciMuSurvive/I");
    tout->Branch("fNeutrinoParentPDG",&fNeutrinoParentPDG,"fNeutrinoParentPDG/I");
    tout->Branch("fNeutrinoParentType",&fNeutrinoParentType,"fNeutrinoParentType/I");
    tout->Branch("fNeutrinoType",&fNeutrinoType,"fNeutrinoType/I");
    tout->Branch("fNeutMode",&fNeutMode,"fNeutMode/I");
    tout->Branch("fNeutrinoParentDecPointSim", "TVector3", &fNeutrinoParentDecPointSim, 128000, spl);
    
    tout->Branch("fMultiplicity",&fMultiplicity,"fMultiplicity/I");
    
    tout->Branch("fMuonNSim",&fMuonNSim,"fMuonNSim/I");
    tout->Branch("fMuonNRec",&fMuonNRec,"fMuonNRec/I");
    tout->Branch("fProtonNSim",&fProtonNSim,"fProtonNSim/I");
    tout->Branch("fProtonNRec",&fProtonNRec,"fProtonNRec/I");
    
    tout->Branch("fMeanBetaRec", &fMeanBetaRec, "fMeanBetaRec/D");
    tout->Branch("fMeanGammaRec", &fMeanGammaRec, "fMeanGammaRec/D");
    
    return tout;
}

TList * NeutrinoUtils::GetHistList()
{
    TList *lout=new TList;
    lout->SetName("lout");
    
    TH1D * hcount = new TH1D("count","",30, -0.5, 29.5); lout->AddAt(hcount,kcount);  
    return lout;
}
